library(tidyverse)
library(sf)
library(terra)
library(ggspatial)

ssp245 = read.csv(file.choose())

ssp245 = ssp245[-c(1,4)]

ssp245$scenario = 'SSP245'

ssp245$period = ifelse(ssp245$period=='2041-2060','2050s',
                       ifelse(ssp245$period=='2061-2080','2070s','Current'))

ssp245$period = factor(ssp245$period,levels=c('Current','2050s','2070s'))

ssp585 = read.csv(file.choose())

ssp585 = ssp585[-c(1,4)]

ssp585$scenario = 'SSP585'

ssp585$period = ifelse(ssp585$period=='2041-2060','2050s',
                       ifelse(ssp585$period=='2061-2080','2070s','Current'))

ssp585$period = factor(ssp585$period,levels=c('Current','2050s','2070s'))

# load shapefile
site1 = read_sf(file.choose()) %>% filter (District_N=='Kyerwa')
# sites are Mbinga and Kyerwa districts


p1 = ggplot()+geom_raster(data=filter(ssp245,GCM=='current'| GCM=='Mean'),aes(x=x,y=y,fill=binary))+
  geom_sf(data=site1,fill='transparent')+
  theme_bw()+
  labs(x='',y='',fill='Suitability classes',title = 'SSP2-4.5')+
  scale_fill_manual(values = c('green4','gray97'))+
  annotation_north_arrow(style=north_arrow_minimal,location='tl',height = unit(0.5,'cm'),width= unit(0.5,'cm'))+
  annotation_scale(location='bl')+# for Arabica (br)
  theme(legend.position = 'none')+
  #facet_wrap(~period,ncol=4)+scale_x_continuous(breaks = seq(34.7,353.0,0.4))# only for Arabica
  facet_wrap(~period,ncol=4)+
  scale_x_continuous(breaks = seq(30.5,31.0,0.2))+# for Robusta
  scale_y_continuous(limits = c(-1.7,-0.9))



p2 = ggplot()+geom_tile(data=filter(ssp585,GCM=='current'| GCM=='Mean'),aes(x=x,y=y,fill=binary))+
  geom_sf(data=site1,fill='transparent')+
  theme_bw()+
  labs(x='',y='',fill='Coffee suitability classes',title = 'SSP5-8.5')+
  scale_fill_manual(values = c('green4','gray97'))+
  annotation_north_arrow(style=north_arrow_minimal,location='tl',height = unit(0.5,'cm'),width= unit(0.5,'cm'))+
  annotation_scale(location='bl')+# for Arabica (br)
  theme(legend.key = element_rect(colour='black',linewidth = 0.75),legend.key.spacing.x = unit(0.03,'in'),legend.position = 'bottom')+
  #facet_wrap(~period,ncol=4)+scale_x_continuous(breaks = seq(34.7,353.0,0.4))+# only for Arabica
  facet_wrap(~period,ncol=4)+
  scale_x_continuous(breaks = seq(30.5,31.0,0.2))+# for Robusta
  guides(fill=guide_legend(title.position='top'))+
  scale_y_continuous(limits = c(-1.7,-0.9)) # for Robusta

p3 = p1/p2

p3

setwd('C:/Users/Mr. Kasongi/Documents/SDM paper revisions and new results')

ggsave('Robusta_ssp245_ssp585_projection.png',p3,height=7.5,width=6,units='in')
#height=8,width=6,units='in'- Arabica

# CONFUSION MATRIX ON THE  SUITABILITY CHANGES
# Arabica nrow=5556
# Arabica nrow=11792
# ssp245
# 2041-2060
mat = ssp245 %>% filter(GCM=='current' | GCM =='Mean' & period =='2050s')

mat1 = data.frame(x=mat$x[1:5556],y=mat$y[1:5556],mat[mat$GCM=='current',]['binary'],mat[mat$GCM=='Mean',]['binary'])

names(mat1)= c('x','y','reference','model')

mat1$reference = factor(mat1$reference,levels=c('Suitable','Unsuitable'))

mat1$model = factor(mat1$model,levels=c('Suitable','Unsuitable'))

# 2061-2080
mat = ssp245 %>% filter(GCM=='current' | GCM =='Mean' & period =='2070s')

mat3 = data.frame(x=mat$x[1:5556],y=mat$y[1:5556],mat[mat$GCM=='current',]['binary'],mat[mat$GCM=='Mean',]['binary'])

names(mat3)= c('x','y','reference','model')

mat3$reference = factor(mat3$reference,levels=c('Suitable','Unsuitable'))

mat3$model = factor(mat3$model,levels=c('Suitable','Unsuitable'))


#------------------------------- NICHE CHANGE-----------------------------------
Future_proj1 = mat1 %>% mutate(category=ifelse(reference=='Suitable' & model=='Suitable',
                                               'Stable',ifelse(reference=='Suitable' & model=='Unsuitable',
                                                'Loss',ifelse(reference=='Unsuitable' & model=='Suitable',
                                                'Gain','none'))),period='2050s')

Future_proj2 = mat3 %>% mutate(category=ifelse(reference=='Suitable' & model=='Suitable',
                                               'Stable',ifelse(reference=='Suitable' & model=='Unsuitable',
                                                'Loss',ifelse(reference=='Unsuitable' & model=='Suitable',
                                                'Gain','none'))),period='2070s')

# merge the dataframe
future_ssp245 = rbind(Future_proj1,Future_proj2)

#-------------------------------------------------------------------------------

# ssp585
# 2041-2060
mat = ssp585 %>% filter(GCM=='current' | GCM =='Mean' & period =='2050s')

mat1 = data.frame(x=mat$x[1:5556],y=mat$y[1:5556],mat[mat$GCM=='current',]['binary'],mat[mat$GCM=='Mean',]['binary'])

names(mat1)= c('x','y','reference','model')

mat1$reference = factor(mat1$reference,levels=c('Suitable','Unsuitable'))

mat1$model = factor(mat1$model,levels=c('Suitable','Unsuitable'))


# 2061-2080
mat = ssp585 %>% filter(GCM=='current' | GCM =='Mean' & period =='2070s')

mat3 = data.frame(x=mat$x[1:5556],y=mat$y[1:5556],mat[mat$GCM=='current',]['binary'],mat[mat$GCM=='Mean',]['binary'])

names(mat3)= c('x','y','reference','model')

mat3$reference = factor(mat3$reference,levels=c('Suitable','Unsuitable'))

mat3$model = factor(mat3$model,levels=c('Suitable','Unsuitable'))


#------------------------------- NICHE CHANGE-----------------------------------
Future_proj1 = mat1 %>% mutate(category=ifelse(reference=='Suitable' & model=='Suitable',
                                               'Stable',ifelse(reference=='Suitable' & model=='Unsuitable',
                                               'Loss',ifelse(reference=='Unsuitable' & model=='Suitable',
                                               'Gain','none'))),period='2050s')

Future_proj2 = mat3 %>% mutate(category=ifelse(reference=='Suitable' & model=='Suitable',
                                               'Stable',ifelse(reference=='Suitable' & model=='Unsuitable',
                                               'Loss',ifelse(reference=='Unsuitable' & model=='Suitable',
                                               'Gain','none'))),period='2070s')

# merge the dataframe
future_ssp585 = rbind(Future_proj1,Future_proj2)


# plot niche change
p5 = ggplot()+geom_raster(data=filter(future_ssp245,category!='none'),aes(x=x,y=y,fill=category))+
  geom_sf(data=site1,fill='transparent')+
  labs(x='',y='',fill='Change in suitability area',title = 'SSP2-4.5')+
  scale_fill_manual(values=c('Loss'='red','Gain'='blue','Stable'='green'))+
  annotation_north_arrow(style=north_arrow_minimal,location='tl',height = unit(0.7,'cm'),width= unit(1,'cm'))+
  annotation_scale(location='br')+
  facet_wrap(~period)+theme_bw()+
  #scale_x_continuous(breaks = seq(34.7,353.0,0.4))+# for Arabica
  facet_wrap(~period,ncol=4)+scale_x_continuous(breaks = seq(30.5,31.0,0.2))+# for Robusta
  theme(legend.position = 'none')+
  scale_y_continuous(limits = c(-1.7,-0.9)) # for Robusta

p6 = ggplot()+geom_raster(data=filter(future_ssp585,category!='none'),aes(x=x,y=y,fill=category))+
  geom_sf(data=site1,fill='transparent')+
  labs(x='',y='',fill='Change in suitability area',title = 'SSP5-8.5')+
  scale_fill_manual(values=c('Loss'='red','Gain'='blue','Stable'='green'))+
  annotation_north_arrow(style=north_arrow_minimal,location='tl',height = unit(0.7,'cm'),width= unit(1,'cm'))+
  annotation_scale(location='br')+
  facet_wrap(~period)+theme_bw()+
  #scale_x_continuous(breaks = seq(34.7,353.0,0.4))+#for Arabica
  facet_wrap(~period,ncol=4)+scale_x_continuous(breaks = seq(30.5,31.0,0.2))+# for Robusta
  theme(legend.position = 'bottom')+
  guides(fill=guide_legend(title.position='top'))+
  scale_y_continuous(limits = c(-1.7,-0.9)) # for Robusta

p7 = p5/p6

p7

# SAVING THE MAP
ggsave('Robusta_niche_change_ssp245_ssp585.png',p7,height=8,width=5,units='in')

############################## THE END #########################################